## pval_cutoff: 0.05
## lfc_cutoff: 1
## low_counts_cutoff: 10

General statistics

# Number of samples
length(counts_data)
## [1] 6
# Number of genes
nrow(counts_data)
## [1] 49315
# Total counts
colSums(counts_data)
## SRR13535276 SRR13535278 SRR13535280 SRR13535288 SRR13535290 SRR13535292 
##     5442647     5284506     6692745    11361755    11378178     4974222

Create DDS objects

# Create DESeqDataSet object
dds <- get_DESeqDataSet_obj(counts_data, ~ treatment)
## [1] TRUE
## [1] TRUE
## [1] "DESeqDataSet object of length 49315 with 0 metadata columns"
## [1] "DESeqDataSet object of length 14816 with 0 metadata columns"
colData(dds)
## DataFrame with 6 rows and 25 columns
##              Assay Type AvgSpotLen       Bases  BioProject    BioSample      Bytes Center Name     Consent DATASTORE filetype DATASTORE provider       DATASTORE region  Experiment treatment GEO_Accession (exp)          Instrument LibraryLayout LibrarySelection  LibrarySource     Organism    Platform                    label ReleaseDate Sample Name            source_name   SRA Study
##             <character>  <numeric>   <numeric> <character>  <character>  <numeric> <character> <character>        <character>        <character>            <character> <character>  <factor>         <character>         <character>   <character>      <character>    <character>  <character> <character>                 <factor>   <POSIXct> <character>            <character> <character>
## SRR13535276     RNA-Seq        300  8225466000 PRJNA694971 SAMN17588686 3252113587         GEO      public          fastq,sra         gs,ncbi,s3 gs.US,ncbi.public,s3..  SRX9943360         A          GSM5043430 Illumina HiSeq 2500        PAIRED             cDNA TRANSCRIPTOMIC Mus musculus    ILLUMINA in space without gravity  2021-09-09  GSM5043430 C2C12 proliferating ..   SRP303354
## SRR13535278     RNA-Seq        300  9203426700 PRJNA694971 SAMN17588684 3619152333         GEO      public          fastq,sra         gs,ncbi,s3 gs.US,ncbi.public,s3..  SRX9943362         A          GSM5043433 Illumina HiSeq 2500        PAIRED             cDNA TRANSCRIPTOMIC Mus musculus    ILLUMINA in space without gravity  2021-09-09  GSM5043433 C2C12 proliferating ..   SRP303354
## SRR13535280     RNA-Seq        300  9323939700 PRJNA694971 SAMN17588682 3735905901         GEO      public          fastq,sra         gs,ncbi,s3 gs.US,ncbi.public,s3..  SRX9943364         A          GSM5043436 Illumina HiSeq 2500        PAIRED             cDNA TRANSCRIPTOMIC Mus musculus    ILLUMINA in space without gravity  2021-09-09  GSM5043436 C2C12 proliferating ..   SRP303354
## SRR13535288     RNA-Seq        300 12863728500 PRJNA694971 SAMN17587373 5128876770         GEO      public          fastq,sra         gs,ncbi,s3 gs.US,ncbi.public,s3..  SRX9943372         C          GSM5043450 Illumina HiSeq 2500        PAIRED             cDNA TRANSCRIPTOMIC Mus musculus    ILLUMINA in space with gravity     2021-09-09  GSM5043450 C2C12 proliferating ..   SRP303354
## SRR13535290     RNA-Seq        300 12849825300 PRJNA694971 SAMN17587371 5136077921         GEO      public          fastq,sra         gs,ncbi,s3 gs.US,ncbi.public,s3..  SRX9943374         C          GSM5043454 Illumina HiSeq 2500        PAIRED             cDNA TRANSCRIPTOMIC Mus musculus    ILLUMINA in space with gravity     2021-09-09  GSM5043454 C2C12 proliferating ..   SRP303354
## SRR13535292     RNA-Seq        300 10569142200 PRJNA694971 SAMN17587369 4229018065         GEO      public          fastq,sra         gs,ncbi,s3 gs.US,ncbi.public,s3..  SRX9943376         C          GSM5043457 Illumina HiSeq 2500        PAIRED             cDNA TRANSCRIPTOMIC Mus musculus    ILLUMINA in space with gravity     2021-09-09  GSM5043457 C2C12 proliferating ..   SRP303354

Sample-to-sample comparisons

# Transform data (blinded rlog)
rld <- get_transformed_data(dds)

PCA plot

pca <- rld$pca
pca_df <- cbind(as.data.frame(colData(dds)) %>% rownames_to_column(var = 'name'), pca$x)
summary(pca)
## Importance of components:
##                            PC1     PC2     PC3     PC4    PC5       PC6
## Standard deviation     32.1475 31.2210 27.8763 24.7058 21.724 4.987e-14
## Proportion of Variance  0.2672  0.2520  0.2009  0.1578  0.122 0.000e+00
## Cumulative Proportion   0.2672  0.5192  0.7202  0.8780  1.000 1.000e+00
ggplot(pca_df, aes(x = PC1, y = PC2, color = label)) +
  geom_point() +
  geom_text(aes(label = name), position = position_nudge(y = -2), show.legend = F, size = 3) +
  scale_color_manual(values = colors_default) +
  scale_x_continuous(expand = c(0.2, 0))

Correlation heatmap

pheatmap(
  cor(rld$matrix),
  annotation_col = as.data.frame(colData(dds)) %>% select(label),
  color = brewer.pal(8, 'YlOrRd')
)

Wald test results

# DE analysis using Wald test
dds_full <- DESeq(dds)
colData(dds_full)
## DataFrame with 6 rows and 26 columns
##              Assay Type AvgSpotLen       Bases  BioProject    BioSample      Bytes Center Name     Consent DATASTORE filetype DATASTORE provider       DATASTORE region  Experiment treatment GEO_Accession (exp)          Instrument LibraryLayout LibrarySelection  LibrarySource     Organism    Platform                    label ReleaseDate Sample Name            source_name   SRA Study sizeFactor
##             <character>  <numeric>   <numeric> <character>  <character>  <numeric> <character> <character>        <character>        <character>            <character> <character>  <factor>         <character>         <character>   <character>      <character>    <character>  <character> <character>                 <factor>   <POSIXct> <character>            <character> <character>  <numeric>
## SRR13535276     RNA-Seq        300  8225466000 PRJNA694971 SAMN17588686 3252113587         GEO      public          fastq,sra         gs,ncbi,s3 gs.US,ncbi.public,s3..  SRX9943360         A          GSM5043430 Illumina HiSeq 2500        PAIRED             cDNA TRANSCRIPTOMIC Mus musculus    ILLUMINA in space without gravity  2021-09-09  GSM5043430 C2C12 proliferating ..   SRP303354   0.667323
## SRR13535278     RNA-Seq        300  9203426700 PRJNA694971 SAMN17588684 3619152333         GEO      public          fastq,sra         gs,ncbi,s3 gs.US,ncbi.public,s3..  SRX9943362         A          GSM5043433 Illumina HiSeq 2500        PAIRED             cDNA TRANSCRIPTOMIC Mus musculus    ILLUMINA in space without gravity  2021-09-09  GSM5043433 C2C12 proliferating ..   SRP303354   0.816941
## SRR13535280     RNA-Seq        300  9323939700 PRJNA694971 SAMN17588682 3735905901         GEO      public          fastq,sra         gs,ncbi,s3 gs.US,ncbi.public,s3..  SRX9943364         A          GSM5043436 Illumina HiSeq 2500        PAIRED             cDNA TRANSCRIPTOMIC Mus musculus    ILLUMINA in space without gravity  2021-09-09  GSM5043436 C2C12 proliferating ..   SRP303354   0.752370
## SRR13535288     RNA-Seq        300 12863728500 PRJNA694971 SAMN17587373 5128876770         GEO      public          fastq,sra         gs,ncbi,s3 gs.US,ncbi.public,s3..  SRX9943372         C          GSM5043450 Illumina HiSeq 2500        PAIRED             cDNA TRANSCRIPTOMIC Mus musculus    ILLUMINA in space with gravity     2021-09-09  GSM5043450 C2C12 proliferating ..   SRP303354   2.071709
## SRR13535290     RNA-Seq        300 12849825300 PRJNA694971 SAMN17587371 5136077921         GEO      public          fastq,sra         gs,ncbi,s3 gs.US,ncbi.public,s3..  SRX9943374         C          GSM5043454 Illumina HiSeq 2500        PAIRED             cDNA TRANSCRIPTOMIC Mus musculus    ILLUMINA in space with gravity     2021-09-09  GSM5043454 C2C12 proliferating ..   SRP303354   1.547156
## SRR13535292     RNA-Seq        300 10569142200 PRJNA694971 SAMN17587369 4229018065         GEO      public          fastq,sra         gs,ncbi,s3 gs.US,ncbi.public,s3..  SRX9943376         C          GSM5043457 Illumina HiSeq 2500        PAIRED             cDNA TRANSCRIPTOMIC Mus musculus    ILLUMINA in space with gravity     2021-09-09  GSM5043457 C2C12 proliferating ..   SRP303354   0.793327
# Wald test results
res <- results(
  dds_full,
  contrast = c('treatment', condition, control),
  alpha = pval_cutoff
)
res
## log2 fold change (MLE): treatment A vs C 
## Wald test p-value: treatment A vs C 
## DataFrame with 14816 rows and 6 columns
##                     baseMean log2FoldChange     lfcSE       stat     pvalue      padj
##                    <numeric>      <numeric> <numeric>  <numeric>  <numeric> <numeric>
## ENSMUSG00000098104   6.35502      0.7602418  1.277274  0.5952064   0.551705  0.916505
## ENSMUSG00000103922   3.00821     -0.4331325  1.505083 -0.2877798   0.773515        NA
## ENSMUSG00000033845 127.62317     -0.3888411  0.565087 -0.6881089   0.491384  0.909368
## ENSMUSG00000102275   2.53432      0.0521924  1.489729  0.0350348   0.972052        NA
## ENSMUSG00000025903  96.02443     -0.2286828  0.361327 -0.6328963   0.526801  0.913742
## ...                      ...            ...       ...        ...        ...       ...
## ENSMUSG00000061654  74.15062     -5.8890703  3.776217 -1.5595158         NA        NA
## ENSMUSG00000079834  36.06458     -0.2951959  0.577612 -0.5110627 0.60930718  0.926139
## ENSMUSG00000095041 247.75852     -0.0310833  0.702407 -0.0442526 0.96470305  0.996240
## ENSMUSG00000063897  25.05491      0.1056571  0.547329  0.1930412 0.84692669  0.978145
## ENSMUSG00000095742   5.67615      3.3261117  1.190578  2.7936945 0.00521097  0.233620
mcols(res)
## DataFrame with 6 rows and 2 columns
##                        type            description
##                 <character>            <character>
## baseMean       intermediate mean of normalized c..
## log2FoldChange      results log2 fold change (ML..
## lfcSE               results standard error: trea..
## stat                results Wald statistic: trea..
## pvalue              results Wald test p-value: t..
## padj                results   BH adjusted p-values
summary(res)
## 
## out of 14816 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up)       : 7, 0.047%
## LFC < 0 (down)     : 60, 0.4%
## outliers [1]       : 159, 1.1%
## low counts [2]     : 2586, 17%
## (mean count < 5)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
plotDispEsts(dds_full)

Summary details

# Upregulated genes (LFC > 0)
res_sig_df %>% filter(log2FoldChange > 0)
# Downregulated genes (LFC < 0)
res_sig_df %>% filter(log2FoldChange < 0)
# Outliers (pvalue and padj are NA)
res[which(is.na(res$pvalue)), ]
## log2 fold change (MLE): treatment A vs C 
## Wald test p-value: treatment A vs C 
## DataFrame with 159 rows and 6 columns
##                     baseMean log2FoldChange     lfcSE      stat    pvalue      padj
##                    <numeric>      <numeric> <numeric> <numeric> <numeric> <numeric>
## ENSMUSG00000103509    7.4013       -5.90608   3.35481  -1.76048        NA        NA
## ENSMUSG00000079554   27.6269       -6.82795   3.46170  -1.97243        NA        NA
## ENSMUSG00000085842   28.7558        4.16157   2.17551   1.91291        NA        NA
## ENSMUSG00000103553   11.1104       -5.49768   2.81475  -1.95317        NA        NA
## ENSMUSG00000047496   20.2028       -2.07168   1.27171  -1.62905        NA        NA
## ...                      ...            ...       ...       ...       ...       ...
## ENSMUSG00000036452  108.1271       -4.32528   1.41232  -3.06254        NA        NA
## ENSMUSG00000053846   27.1020       -5.89974   2.22511  -2.65144        NA        NA
## ENSMUSG00000024867   18.2290       -1.55450   1.46155  -1.06359        NA        NA
## ENSMUSG00000048029   22.7289       -6.52465   3.90834  -1.66942        NA        NA
## ENSMUSG00000061654   74.1506       -5.88907   3.77622  -1.55952        NA        NA
# Low counts (only padj is NA)
res[which(is.na(res$padj) & !is.na(res$pvalue)), ]
## log2 fold change (MLE): treatment A vs C 
## Wald test p-value: treatment A vs C 
## DataFrame with 2586 rows and 6 columns
##                     baseMean log2FoldChange     lfcSE       stat    pvalue      padj
##                    <numeric>      <numeric> <numeric>  <numeric> <numeric> <numeric>
## ENSMUSG00000103922   3.00821     -0.4331325   1.50508 -0.2877798  0.773515        NA
## ENSMUSG00000102275   2.53432      0.0521924   1.48973  0.0350348  0.972052        NA
## ENSMUSG00000103280   2.25921     -0.2036441   1.54506 -0.1318034  0.895140        NA
## ENSMUSG00000103355   3.20134     -0.0484038   1.26490 -0.0382670  0.969475        NA
## ENSMUSG00000103845   2.38594     -1.2816787   1.97493 -0.6489755  0.516354        NA
## ...                      ...            ...       ...        ...       ...       ...
## ENSMUSG00000064347   2.33249       2.646153   1.66191   1.592236  0.111332        NA
## ENSMUSG00000064360   3.81827       0.236955   1.17589   0.201512  0.840298        NA
## ENSMUSG00000064364   4.38183      -1.462709   1.21266  -1.206204  0.227739        NA
## ENSMUSG00000064365   2.11846       1.897276   1.77519   1.068774  0.285171        NA
## ENSMUSG00000064366   4.04585      -1.225589   2.08042  -0.589107  0.555790        NA

Shrunken LFC results

plotMA(res)

# Shrunken LFC results
res_shrunken <- lfcShrink(
  dds_full,
  coef = str_c('treatment_', condition, '_vs_', control),
  type = 'apeglm'
)
res_shrunken
## log2 fold change (MAP): treatment A vs C 
## Wald test p-value: treatment A vs C 
## DataFrame with 14816 rows and 5 columns
##                     baseMean log2FoldChange     lfcSE     pvalue      padj
##                    <numeric>      <numeric> <numeric>  <numeric> <numeric>
## ENSMUSG00000098104   6.35502     0.01967164  0.206079   0.551705        NA
## ENSMUSG00000103922   3.00821    -0.00811498  0.205596   0.773515        NA
## ENSMUSG00000033845 127.62317    -0.04685436  0.201406   0.491384  0.910778
## ENSMUSG00000102275   2.53432     0.00104597  0.205357   0.972052        NA
## ENSMUSG00000025903  96.02443    -0.05793424  0.187805   0.526801  0.914774
## ...                      ...            ...       ...        ...       ...
## ENSMUSG00000061654  74.15062    -0.00862303  0.207574         NA        NA
## ENSMUSG00000079834  36.06458    -0.03367781  0.198685 0.60930718  0.925945
## ENSMUSG00000095041 247.75852    -0.00209616  0.198841 0.96470305  0.995670
## ENSMUSG00000063897  25.05491     0.01339576  0.194495 0.84692669  0.977817
## ENSMUSG00000095742   5.67615     2.01391517  1.919874 0.00521097        NA
plotMA(res_shrunken)

mcols(res_shrunken)
## DataFrame with 5 rows and 2 columns
##                        type            description
##                 <character>            <character>
## baseMean       intermediate mean of normalized c..
## log2FoldChange      results log2 fold change (MA..
## lfcSE               results posterior SD: treatm..
## pvalue              results Wald test p-value: t..
## padj                results   BH adjusted p-values
summary(res_shrunken, alpha = pval_cutoff)
## 
## out of 14816 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up)       : 7, 0.047%
## LFC < 0 (down)     : 60, 0.4%
## outliers [1]       : 159, 1.1%
## low counts [2]     : 3159, 21%
## (mean count < 6)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results

Summary details

# Upregulated genes (LFC > 0)
res_shrunken_sig_df %>% filter(log2FoldChange > 0)
# Downregulated genes (LFC < 0)
res_shrunken_sig_df %>% filter(log2FoldChange < 0)
# Outliers (pvalue and padj are NA)
res_shrunken[which(is.na(res_shrunken$pvalue)), ]
## log2 fold change (MAP): treatment A vs C 
## Wald test p-value: treatment A vs C 
## DataFrame with 159 rows and 5 columns
##                     baseMean log2FoldChange     lfcSE    pvalue      padj
##                    <numeric>      <numeric> <numeric> <numeric> <numeric>
## ENSMUSG00000103509    7.4013     -0.0123157  0.207854        NA        NA
## ENSMUSG00000079554   27.6269     -0.0114825  0.207787        NA        NA
## ENSMUSG00000085842   28.7558      0.0246738  0.209305        NA        NA
## ENSMUSG00000103553   11.1104     -0.0177780  0.208382        NA        NA
## ENSMUSG00000047496   20.2028     -0.0492316  0.214096        NA        NA
## ...                      ...            ...       ...       ...       ...
## ENSMUSG00000036452  108.1271    -0.06071289  0.220435        NA        NA
## ENSMUSG00000053846   27.1020    -0.03057686  0.210614        NA        NA
## ENSMUSG00000024867   18.2290    -0.02871114  0.208606        NA        NA
## ENSMUSG00000048029   22.7289    -0.00870315  0.207584        NA        NA
## ENSMUSG00000061654   74.1506    -0.00862303  0.207574        NA        NA
# Low counts (only padj is NA)
res_shrunken[which(is.na(res_shrunken$padj) & !is.na(res_shrunken$pvalue)), ]
## log2 fold change (MAP): treatment A vs C 
## Wald test p-value: treatment A vs C 
## DataFrame with 3159 rows and 5 columns
##                     baseMean log2FoldChange     lfcSE     pvalue      padj
##                    <numeric>      <numeric> <numeric>  <numeric> <numeric>
## ENSMUSG00000098104   6.35502     0.01967164  0.206079   0.551705        NA
## ENSMUSG00000103922   3.00821    -0.00811498  0.205596   0.773515        NA
## ENSMUSG00000102275   2.53432     0.00104597  0.205357   0.972052        NA
## ENSMUSG00000103280   2.25921    -0.00361950  0.205495   0.895140        NA
## ENSMUSG00000103355   3.20134    -0.00120582  0.204580   0.969475        NA
## ...                      ...            ...       ...        ...       ...
## ENSMUSG00000065947   5.53967      0.0697355  0.221195 0.07681665        NA
## ENSMUSG00000064364   4.38183     -0.0419282  0.210623 0.22773891        NA
## ENSMUSG00000064365   2.11846      0.0237202  0.208145 0.28517145        NA
## ENSMUSG00000064366   4.04585     -0.0114088  0.206890 0.55578956        NA
## ENSMUSG00000095742   5.67615      2.0139152  1.919874 0.00521097        NA

Visualizing results

Heatmaps

# Plot normalized counts (z-scores)
pheatmap(counts_sig_norm[2:7], 
         color = brewer.pal(8, 'YlOrRd'), 
         cluster_rows = T, 
         show_rownames = F,
         annotation_col = as.data.frame(colData(dds)) %>% select(label),
         border_color = NA,
         fontsize = 10,
         scale = 'row',
         fontsize_row = 10, 
         height = 20)

# Plot log-transformed counts
pheatmap(counts_sig_log[2:7], 
         color = rev(brewer.pal(8, 'RdYlBu')), 
         cluster_rows = T, 
         show_rownames = F,
         annotation_col = as.data.frame(colData(dds)) %>% select(label),
         border_color = NA,
         fontsize = 10,
         fontsize_row = 10, 
         height = 20)

# Plot log-transformed counts (top 24 DE genes)
pheatmap(counts_sig_log %>% filter(ensembl_gene_id %in% (res_sig_df %>% head(24))$ensembl_gene_id) %>% select(-ensembl_gene_id) %>% column_to_rownames(var = 'mgi_symbol'),
         color = rev(brewer.pal(8, 'RdYlBu')), 
         cluster_rows = T, 
         show_rownames = T,
         annotation_col = as.data.frame(colData(dds)) %>% select(label), 
         fontsize = 10,
         fontsize_row = 10, 
         height = 20)

Volcano plots

# Unshrunken LFC
res_df %>% 
  mutate(
    sig_threshold = if_else(
      padj < pval_cutoff & abs(log2FoldChange) >= lfc_cutoff,
      if_else(log2FoldChange > 0, 'DE-up', 'DE-down'),
      'non-DE'
    )
  ) %>% 
  filter(!is.na(sig_threshold)) %>% 
  ggplot() +
  geom_point(aes(x = log2FoldChange, y = -log10(padj), colour = sig_threshold)) +
  scale_color_manual(values = c('blue', 'red', 'gray')) +
  xlab('log2 fold change') + 
  ylab('-log10 adjusted p-value')

# Shrunken LFC
res_shrunken_df %>% 
  mutate(
    sig_threshold = if_else(
      padj < pval_cutoff & abs(log2FoldChange) >= lfc_cutoff,
      if_else(log2FoldChange > 0, 'DE-up', 'DE-down'),
      'non-DE'
    )
  ) %>% 
  filter(!is.na(sig_threshold)) %>% 
  ggplot() +
  geom_point(aes(x = log2FoldChange, y = -log10(padj), colour = sig_threshold)) +
  scale_color_manual(values = c('blue', 'red', 'gray')) +
  xlab('log2 fold change') + 
  ylab('-log10 adjusted p-value')

GSEA (all)

Hallmark genesets

# Shrunken LFC
get_fgsea_res(rank_lfc, mm_h) %>% plot_enrichment_table(rank_lfc, mm_h)

# Wald stat
get_fgsea_res(rank_stat, mm_h) %>% plot_enrichment_table(rank_stat, mm_h)

# Rank: sign(LFC) * -log10(pvalue)
get_fgsea_res(rank_pval, mm_h) %>% plot_enrichment_table(rank_pval, mm_h)

GO biological process

# Shrunken LFC
get_fgsea_res(rank_lfc, mm_c5_bp) %>% plot_enrichment_table(rank_lfc, mm_c5_bp)

# Wald stat
get_fgsea_res(rank_stat, mm_c5_bp) %>% plot_enrichment_table(rank_stat, mm_c5_bp)

# Rank: sign(LFC) * -log10(pvalue)
get_fgsea_res(rank_pval, mm_c5_bp) %>% plot_enrichment_table(rank_pval, mm_c5_bp)

GO cellular component

# Shrunken LFC
get_fgsea_res(rank_lfc, mm_c5_cc) %>% plot_enrichment_table(rank_lfc, mm_c5_cc)

# Wald stat
get_fgsea_res(rank_stat, mm_c5_cc) %>% plot_enrichment_table(rank_stat, mm_c5_cc)

# Rank: sign(LFC) * -log10(pvalue)
get_fgsea_res(rank_pval, mm_c5_cc) %>% plot_enrichment_table(rank_pval, mm_c5_cc)

GO molecular function

# Shrunken LFC
get_fgsea_res(rank_lfc, mm_c5_mf) %>% plot_enrichment_table(rank_lfc, mm_c5_mf)

# Wald stat
get_fgsea_res(rank_stat, mm_c5_mf) %>% plot_enrichment_table(rank_stat, mm_c5_mf)

# Rank: sign(LFC) * -log10(pvalue)
get_fgsea_res(rank_pval, mm_c5_mf) %>% plot_enrichment_table(rank_pval, mm_c5_mf)

GSEA (DE)

Hallmark genesets

# Shrunken LFC
get_fgsea_res(rank_lfc, mm_h) %>% plot_enrichment_table(rank_lfc, mm_h)

# Wald stat
get_fgsea_res(rank_stat, mm_h) %>% plot_enrichment_table(rank_stat, mm_h)

# Rank: sign(LFC) * -log10(pvalue)
get_fgsea_res(rank_pval, mm_h) %>% plot_enrichment_table(rank_pval, mm_h)

GO biological process

# Shrunken LFC
get_fgsea_res(rank_lfc, mm_c5_bp) %>% plot_enrichment_table(rank_lfc, mm_c5_bp)

# Wald stat
get_fgsea_res(rank_stat, mm_c5_bp) %>% plot_enrichment_table(rank_stat, mm_c5_bp)

# Rank: sign(LFC) * -log10(pvalue)
get_fgsea_res(rank_pval, mm_c5_bp) %>% plot_enrichment_table(rank_pval, mm_c5_bp)

GO cellular component

# Shrunken LFC
get_fgsea_res(rank_lfc, mm_c5_cc) %>% plot_enrichment_table(rank_lfc, mm_c5_cc)

# Wald stat
get_fgsea_res(rank_stat, mm_c5_cc) %>% plot_enrichment_table(rank_stat, mm_c5_cc)

# Rank: sign(LFC) * -log10(pvalue)
get_fgsea_res(rank_pval, mm_c5_cc) %>% plot_enrichment_table(rank_pval, mm_c5_cc)

GO molecular function

# Shrunken LFC
get_fgsea_res(rank_lfc, mm_c5_mf) %>% plot_enrichment_table(rank_lfc, mm_c5_mf)

# Wald stat
get_fgsea_res(rank_stat, mm_c5_mf) %>% plot_enrichment_table(rank_stat, mm_c5_mf)

# Rank: sign(LFC) * -log10(pvalue)
get_fgsea_res(rank_pval, mm_c5_mf) %>% plot_enrichment_table(rank_pval, mm_c5_mf)

System info

sessionInfo()
## R version 4.1.0 (2021-05-18)
## Platform: x86_64-conda-linux-gnu (64-bit)
## Running under: CentOS Linux 7 (Core)
## 
## Matrix products: default
## BLAS/LAPACK: /home/chan/mRNA_seq_pipeline/.snakemake/conda/9a19315a020c824d12f8055f7c009b0f/lib/libopenblasp-r0.3.18.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8     LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8    LC_PAPER=en_US.UTF-8       LC_NAME=C                  LC_ADDRESS=C               LC_TELEPHONE=C             LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] fgsea_1.20.0                RColorBrewer_1.1-2          pheatmap_1.0.12             DESeq2_1.34.0               SummarizedExperiment_1.24.0 Biobase_2.54.0              MatrixGenerics_1.6.0        matrixStats_0.61.0          GenomicRanges_1.46.0        GenomeInfoDb_1.30.0         IRanges_2.28.0              S4Vectors_0.32.0            BiocGenerics_0.40.0         scales_1.1.1                forcats_0.5.1               stringr_1.4.0               dplyr_1.0.7                 purrr_0.3.4                 readr_2.1.1                 tidyr_1.1.4                 tibble_3.1.6                ggplot2_3.3.5               tidyverse_1.3.1            
## 
## loaded via a namespace (and not attached):
##  [1] colorspace_2.0-2       ellipsis_0.3.2         XVector_0.34.0         fs_1.5.1               rstudioapi_0.13        farver_2.1.0           bit64_4.0.5            mvtnorm_1.1-3          AnnotationDbi_1.56.1   fansi_0.4.2            apeglm_1.16.0          lubridate_1.8.0        xml2_1.3.3             splines_4.1.0          cachem_1.0.6           geneplotter_1.72.0     knitr_1.35             jsonlite_1.7.2         broom_0.7.10           annotate_1.72.0        dbplyr_2.1.1           png_0.1-7              compiler_4.1.0         httr_1.4.2             backports_1.4.0        assertthat_0.2.1       Matrix_1.3-4           fastmap_1.1.0          cli_3.1.0              htmltools_0.5.2        tools_4.1.0            coda_0.19-4            gtable_0.3.0           glue_1.5.1             GenomeInfoDbData_1.2.7 fastmatch_1.1-3        Rcpp_1.0.7             bbmle_1.0.24           cellranger_1.1.0       jquerylib_0.1.4        vctrs_0.3.8            Biostrings_2.62.0      xfun_0.28              rvest_1.0.2            lifecycle_1.0.1        XML_3.99-0.8           MASS_7.3-54            zlibbioc_1.40.0        vroom_1.5.7            hms_1.1.1              parallel_4.1.0         yaml_2.2.1             memoise_2.0.1          gridExtra_2.3          emdbook_1.3.12         bdsmatrix_1.3-4        stringi_1.7.6          RSQLite_2.2.8          highr_0.9              genefilter_1.76.0      BiocParallel_1.28.0    rlang_0.4.12           pkgconfig_2.0.3        bitops_1.0-7           evaluate_0.14          lattice_0.20-45        labeling_0.4.2         bit_4.0.4              tidyselect_1.1.1       plyr_1.8.6             magrittr_2.0.1         R6_2.5.1               generics_0.1.1         DelayedArray_0.20.0    DBI_1.1.1              pillar_1.6.4           haven_2.4.3            withr_2.4.3            survival_3.2-13        KEGGREST_1.34.0        RCurl_1.98-1.5         modelr_0.1.8           crayon_1.4.2           utf8_1.2.2             tzdb_0.2.0             rmarkdown_2.11         locfit_1.5-9.4         grid_4.1.0             readxl_1.3.1           data.table_1.14.2      blob_1.2.2             reprex_2.0.1           digest_0.6.29          xtable_1.8-4           numDeriv_2016.8-1.1    munsell_0.5.0